Variational inference for Markov jump processes

Markov jump processes play an important role in a large number of application domains. However, realistic systems are analytically intractable and they have traditionally been analysed using simulation based techniques, which do not provide a framework for statistical inference. We propose a mean field approximation to perform posterior inference and parameter estimation. The approximation allows a practical solution to the inference problem, while still retaining a good degree of accuracy. We illustrate our approach on two biologically motivated systems. Introduction Markov jump processes (MJPs) underpin our understanding of many important systems in science and technology. They provide a rigorous probabilistic framework to model the joint dynamics of groups (species) of interacting individuals, with applications ranging from information packets in a telecommunications network to epidemiology and population levels in the environment. These processes are usually non-linear and highly coupled, giving rise to non-trivial steady states (often referred to as emerging properties). Unfortunately, this also means that exact statistical inference is unfeasible and approximations must be made in the analysis of these systems. A traditional approach, which has been very successful throughout the past century, is to ignore the discrete nature of the processes and to approximate the stochastic process with a deterministic process whose behaviour is described by a system of non-linear, coupled ODEs. This approximation relies on the stochastic fluctuations being negligible compared to the average population counts. There are many important situations where this assumption is untenable: for example, stochastic fluctuations are reputed to be responsible for a number of important biological phenomena, from cell differentiation to pathogen virulence [1]. Researchers are now able to obtain accurate estimates of the number of macromolecules of a certain species within a cell [2, 3], prompting a need for practical statistical tools to handle discrete data. Sampling approaches have been extensively used to simulate the behaviour of MJPs. Gillespie’s algorithm and its generalisations [4, 5] form the basis of many simulators used in systems biology studies. The simulations can be viewed as individual samples taken from a completely specified MJP, and can be very useful to reveal possible steady states. However, it is not clear how observed data can be incorporated in a principled way, which renders this approach of limited use for posterior inference and parameter estimation. A Markov chain Monte Carlo (MCMC) approach to incorporate observations has been recently proposed by Boys et al. [6]. While this approach holds a lot of promise, it is computationally very intensive. Despite several simplifying approximations, the correlations between samples mean that several millions of MCMC iterations are needed even in simple examples. In this paper we present an alternative, deterministic approach to posterior inference and parameter estimation in MJPs. We extend the mean-field (MF) variational approach ([cf. e.g. 7]) to approximate a probability distribution over an (infinite dimensional) space of discrete paths, representing the time-evolving state of the system. In this way, we replace the couplings between the 1 different species by their average, mean-field (MF) effect. The result is an iterative algorithm that allows parameter estimation and prediction with reasonable accuracy and very contained computational costs. The rest of this paper is organised as follows: in sections 1 and 2 we review the theory of Markov jump processes and introduce our general strategy to obtain a MF approximation. In section 3 we introduce the Lotka-Volterra model which we use as an example to describe how our approach works. In section 4 we present experimental results on simulated data from the Lotka-Volterra model and from a simple gene regulatory network. Finally, we discuss the relationship of our study to other stochastic models, as well as further extensions and developments of our approach. 1 Markov jump processes We start off by establishing some notation and basic definitions. A D-dimensional discrete stochastic process is a family of D-dimensional discrete random variables x(t) indexed by the continuous time t. In our examples, the values taken by x(t) will be restricted to the non-negative integers N D 0 . The dimensionality D represents the number of (molecular) species present in the system; the components of the vector x (t) then represent the number of individuals of each species present at time t. Furthermore, the stochastic processes we will consider will always be Markovian, i.e. given any sequence of observations for the state of the system (xt1 , . . . ,xtN ), the conditional probability of the state of the system at a subsequent time xtN+1 depends only on the last of the previous observations. A discrete stochastic process which exhibits the Markov property is called a Markov jump process (MJP). A MJP is characterised by its process rates f (x|x), defined ∀x 6= x; in an infinitesimal time interval δt, the quantity f (x|x) δt represents the infinitesimal probability that the system will make a transition from state x at time t to state x at time t + δt. Explicitly, p (x|x) ' δx′x + δtf (x |x) (1) where δx′x is the Kronecker delta and the equation becomes exact in the limit δt → 0. Equation (1) implies by normalisation that f (x|x) = − ∑ x ′ 6=x f (x |x). The interpretation of the process rates as infinitesimal transition probabilities highlights the simple relationship between the marginal distribution pt (x) and the process rates. The probability of finding the system in state x at time t + δt will be given by the probability that the system was already in state x at time t, minus the probability that the system was in state x at time t and jumped to state x, plus the probability that the system was in a different state x′′ at time t and then jumped to state x. In formulae, this is given by pt+δt (x) = pt (x) 